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ABSTRACT 


In this paper, we present closed-form formulas for the lift and moment coefficients of a lifting surface 
in two-dimensional, unsteady, compressible, subsonic flow utilizing a newly developed explicit 
analytical solution of the Possio equation. Numerical calculations are consistent with previous numerical 
tables based on series expansions or ad hoc numerical schemes. More importantly, these formulas lend 
themselves readily to flutter analysis, compared with the tedious table-look-up schemes currently in use. 
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CT 


dummy variable 
<|> velocity potential 

X dummy variable 

to frequency of oscillation, rad/sec 


1. INTRODUCTION 

The pressure distribution on a lifting surface in unsteady aerodynamics and the lift/moment 
coefficients deduced therefrom are essential for many aerodynamic/aeroelastic calculations including 
flutter analysis. A central role is played in this theory by the Possio (integral) equation, relating the 
normal velocity of the fluid (downwash) to the pressure distribution of a lifting surface in two- 
dimensional, oscillatory, subsonic compressible flow, derived by Possio. 1 Despite the effort of many 
aerodynamicists and aeroelasticians, however, no explicit solution has so far been found. Dietze and 
later Fettis 3 and others turned to creating tables of lift and moment coefficients based on various 
numerical approximations. The lack of explicit formulas for the aerodynamic coefficients, in turn, 
resulted in making the aeroelastic calculations a tedious iterative process. Haskind, 4 Reissner, 5 and 
Timman, et al., 6 formulated the two-dimensional flow problem in elliptic coordinates and this resulted in 
a series expansion in terms of Mathieu functions and also developed aerodynamic coefficient tables. By 
the 1950’s, the aerodynamic theory had extended to the more general three-dimensional lifting surface 
problem. Watkins, Runyan and Woolston 7 used series expansions which led to a numerical scheme 
known as the Kernel Method (KFM). In the late 1960’s, Albano and Rodden 8 developed a numerical 
scheme called the “doublet-lattice method” which extends to even more complicated nonplanar 
configurations, but without theoretical justification. 8, 9 As the flow problem became more complicated it 
became even more difficult to combine the structure with the aerodynamics necessary for the aeroelastic 
problem. 

Recently, Balakrishnan 10 derived an explicit solution to the Possio equation, which is accurate to the 

2 

order of M 7 iog M . This approximation is certainly valid for M<0.7, above which the linearized 
l-M“ l-M 2 2 

theory is generally conceded as no longer valid. Figure 1 shows that M 7 iog — is well below one for 

I - M“ l-M 2 

M < 0.7 . In this paper we derive explicit formulas for the lift and moment for any two-dimensional 
lifting surface in subsonic compressible flow, based on the Balakrishnan solution. 


We begin in section 2 with the Laplace transform version of the Possio equation and its solution 
derived by Balakrishnan. 10 In section 3, we proceed to obtain closed-form formulas for the lift and 
moment. Numerical calculations are given in section 4, which are compared with results obtained by 

*1 (L \ \ 

previous authors, ' ’ ’ showing acceptable agreement. In section 5 we illustrate the ease with which 
the lift and moment formulas are applied to flutter analysis as compared with table-look-up schemes, by 
considering a standard model for wings with large aspect ratio. Concluding remarks are in the final 
section. 
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Figure 1. Plot of -!^log— and M 2 logM. 
1 - M 1 - M 


2. SOLUTION OF POSSIO EQUATION - LAPLACE TRANSFORM VERSION 

In this section we present the Laplace transform version of the Possio equation and its solution, 
referring to reference 10 for more detailed analysis. To clarify the terminology, we begin with the two- 
dimensional field equations for the velocity potential in (linearized) compressive flow: 

1 a 1 2 <l>(jc, z, r) t 2 L , a 2 (t>(*,z,r) p2 a 2 <l>(jc, z, r) 3 2 <K*, z, t) _ 0 (2.1) 

a 2 d 2 dx 2 dz 2 

“oo 

-oo < x < °°, 0 < Z < 


subject to the boundary conditions: 

i) Flow Tangency 


3(|)(x, 0 + , 0 

dz 


w a (x, t), |x| < b\ 


( 2 . 2 ) 


ii) Kutta-Joukowski 


3<()(x, 0 , t) + ^ d<j>(x, 0 , r) _ q, x = ft - and x > b', (2.3) 

dt dx 

iii) and the usual vanishing far-field conditions. 
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A “particular” solution of (2.1), which satisfies iii) and ii) for x> b is given as (see ref. 12) 




1 + 


x jc-C 5L + P 2 (n 2 ± * 2 ) 


t/p 


2 (/ 


«~P 


dX. (2-4) 


/~2 D 2 , 2 2 T 

VX + P (T| +Z ) 

It is a “particular” solution in that the initial conditions are assumed to be: 

x/ ™ _ 3<K*, z.0) _ n 

<|)(x,z,0) = ^ U, 

which are necessary for the Laplace transform theory. To satisfy the flow-tangency condition, we 
substitute this solution into (2.2), and must obtain 


w (x, t) = li m -r - f < ^*1 J 

flV 47t J -b J -~ J -~ d z 2 


x *-C 7x 2 + pW + zV 


t*+ 


t/p 


2 G 


«-P 


Vx 2 +pW+z 2 ) 


rfx. (2.5) 


Next, we take the Laplace transform on both sides by defining 

w a (x, X) = J” e~ X ‘w a (x, t)dt ; A(x, X) = j^ e A(x,t)dt, ReX>0. 


We now have: 




exp 


f/P" 

f Vx 2 P 2 (~n 2 + ^ 2 ) 


3" r 00 

3z 2 


«-P 


#2 n2) 2 2. 

VX + P Ol +2 ) 




(2.6) 


U p 


2 (/ 


A(C,X)- lim — ^A'qI 

x-^o 3z 


7x 2 + p 2 ^ 


«-P 


This leads to (see ref. 10) the “Laplace transform version” of the Possio equation, 

w(x,X) = —X (G(jc-C,M,X)AP(C,X)JC, 
pt /' J ~ b 


(2.7) 
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where 


A P(C,X) = -p (M(C,X), 



G(y, M, X) is analytic in X for the entire complex plane except for the “branch cut” along the negative 
real axis, due to the essential singularity of the modified Bessel functions. Hence, we can set X = id ) , 
yielding the usual Possio equation for “oscillatory motion.” 13 

The explicit solution to (2.7) developed by Balakrishnan in ref. 10 is given by 

f(x,X) = r(j:, X) + coshy(x- a)r(a, X)dcr + sinhY(x:- o)r(CT, X)rfo (2.8) 


where 


f(x, X) = g(x, X)-yK(y,x)j ^ sinhy(l -a)r( a, X)do 

+ (£p(p, x)-MyX"(y, jc))J coshY(l -a)r(o, X)da, 


(2.8a) 


and 


f(x,X) 
g(x, X) 

K(y,x) 

p(^x) 


e Myx AP(x, X), 

2 / 1 — jr r 1 + r e My WX) 

c-, * 


i fiTc W-0) 


i/Eif' II 

n 2 y 1 +xJ-i V 1 


C-JC 




ll-Xl 

r,-v 1 j 

f2 + i 

h+xJ 

'o x-t- 1 

t 


dt. 


For ease of notation, we have set b = 1 since it is only a matter of scaling. The coefficients 


r 1 r 1 

J sinhy( 1 - <J)r(o)do and J coshy( 1 - o)r(o)do 
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are readily calculated. Let 


Sh(q) = J 1 ^ Y sinhy( 1 -o)q(a)do 
Ch(q ) = J 1 coshy(l -o)q(o)da. 


Then from (2.8a) we obtain the linear equations: 

Sh(r) = Sh(g ) - Sh(K)Sh(r) + [kSh(p) - UySh(K)]Ch(r) ( 2 . 9 ) 

Ch(r) = Ch(g) - Ch(K)Sh(r) + [kCh(p) - MyCh(K)]Ch(r). 


The determinant is 


D = (1 + Sh(K))(l -kCh(p)) + (My +kSh(p))Ch(K). 


A plot of D as a function of the reduced frequency is given in figure 2 for a family of values of M < 0.7 
showing it is actually positive for these values, so that (2.9) has a unique solution. 



Figure 2. The determinant D as a function of the reduced frequency (k) for various Mach numbers. 


For the incompressible case, the Balaknshnan solution simplifies this and we obtain 

f(x,X) = r(x, X) + r(a, X)da 


( 2 . 10 ) 
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and 


where 


r (*) = g(*) + kp(k, .v)J ^ r(o)da , 


ri I , 8 ^ da i -* 1 

J r(a)da = = = J- g(a)da 

J -' 1 -rf p(k,a)do kK 0 (k) + K 


( 2 . 11 ) 


1 r 1 


D(k ) J - 1 


J g{o)da. 


Plugging (2.11) into (2.10), we have 


f(x,X) = g(*) + &r g(o)da + k — — f* g{o)da\ p(k, x) + kf p(k,a)do 
-1 D(k) J -i L J -l 

which agrees with the known (e.g. Kiissner-Schwarz 14 ) solution for “oscillatory motion. 


( 2 . 12 ) 


3. LIFT AND MOMENT CALCULATIONS 

In this section we calculate the lift and moment based on Balakrishnan’s solution, equation (2.8), 


accurate again to the order of -log — 52— (or equivalently M^logM ). 

1-M 1-M~ 


The lift is given by: 


L(X) = -f 1 A P{x, X)dx 

J -l 

= -J_, ''(*) + ltJ ( coshy(j:-a)r(a)(ia + Yj a sinh y(x-o)r(cf)da 


dx 


(3.1) 


My 

■ £ » (r) - e M7 M(r) - 


where 


E 0 (r) = f‘ e Myx r(x)dx. 
« _ 1 
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Similarly we have for the moment about the axis x - a as 


M(k) = ^(x-a)AP(x,\)dx 

= J 1 e MYJt (x-a)|V(.x) + pJ^coshY(*-oM°M CT 

+ yJ* sinhY(x-a)r(a)daJc?j: (3.2) 

= f-^- a ) £ o <r) + £|(r) '^i eMTc ' l( '' ) 

k J Mi 

+ -i-fl + i-ale MY5 ^ r >’ 

Myv k ) 

where 

£,(r) = j' ^ xe Myx r(x)dx. 

Hence to calculate lift and moment we need to calculate two more constants namely, 

E 0 (r) and E,(r). 

We have again from (2.8a), 

E 0 (r) = E 0 (g)-E 0 (K)Sh(r)+[kE 0 (p)-MyE 0 (K)]Ch(r), 

E x {r) = E x (g) - E x (K)Sh(r) + [kE x (p) -MyE x (K)]Ch(r). 

We are then left to find the twelve constants, E 0 (g), E x (g ), Sh(g), Ch(g), E 0 (p), E x (p), Sh(p ), 
Ch(p ) , E 0 (K) ,£,(£), Sh(K) and Ch{K). It is difficult to obtain explicit formulas for these constants. 

But since we are only interested in those terms of the order M 2 logM , we expand them in a power series 

of M. Keeping terms up to M logM , we obtain 

L(k) = L 0 (^) + (M 2 logM)L M (?i) and M (A.) = M 0 (X) + (M 2 logM)M M (X), 
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where 


L 0 (A) = ~2pUj^ J^[k(l-o) + C(k)]w a (o,k)do, 

M 0 (X) = p(/J - 2ak- (1 + 2a)C{k) + [2 + *(1 + 2a)]a -ko 2 }\v a (a, A)jdo, 

£m(*) = -P^J 1 ,^^^ 3 +3k 2 C(k) + 2kC(k) 2 - a[& 3 + 2 k 2 C(fc)]}w a (c, X)do, (33) 

M m (X) = p(/f' /i±2{- a it 3 -(l +3a) /t 2 C(/t)-(l +2a) kC(k) 2 
j_i a/ i — a 

+ a[ak 3 + (1 + 2a) k 2 C{k)] }w fl (a, X)do 

and 


C(k) = 


K ,(t) 

*,(*) + *<,(*) 


Theodorsen function. 


Equation (3.3) is for a general class of lifting surfaces, in that the downwash function is arbitrary. For 
a typical section we can express the downwash as 


w a (x,t) = -h(t) -(x-a)a(t)-Ua.(t), for \x\<b, 
or 

w a (A) = -h(X) - (x - a)Xa(X) - Ua(X). (3.4) 


We ignore the initial conditions here since we will be interested only in the stability of the wing. We 
could also easily include control surfaces, but for simplicity we will not do so. 

Plugging (3.4) into (3.3) we obtain: 


L(X) = npU 2 b[L k Ck)^ + L a Ck)a(X)l 
M(X) = npU 2 b 2 [M h Ck)^ + M a Ck)a(X)\ 


(3.5) 
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where 


LJk) 


k 2 + 2kC(k) + M 2 logM 


k 4 


2k 3 C(k) + 2k 2 C(ic) 2 ^, 


k-ak 2 + C(k)[2 + (l-2a)k] 
fr3 H 

+ M 2 logMj — — a + C(k ) 


2k 2 + (^~ 2ajk 3 + C(k) 2 [2k + ( 1 - 2a)k 2 ] 


M h (k) = ait 2 + (1 + 2a)kC(k) + M 2 logMi (1 + 2a)k 2 C{k) 2 + Q + 2a)* 3 C{k) + ^-a L 


M a (k) = [a - A)* + a 2 )* 2 + C(*)[l +2 a + Q - 2 a 2 )*] + M 2 logM ?* 3 - a -k 4 


+ C(k) 2 


(1 +2a)* + | *-2a 2 )* 2 + C(*) (^ + 2a)* 2 -2a 2 * 3 


Smilg and Wassermann 15 defined the aerodynamic coefficients differently, they had 


L(X) = + 


M(k) = -npU 2 b 2 k 2 \ 




h(X) 


~M a k M + a)(Z a (*) + M h (k )) + Q + a)\(*)]a(X) 


(3.6) 
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Here the coefficients are given by: 


L h (k ) = 1 + lc(k) + M 2 logM 

k 


k 2 - - - 21 

y + 2kC(k) + 2C(k) I, 


L a (k ) = \ + l + C(k) 
k 


2 _ 2 
l 2 


k 2 k. 


+ M 2 logM| \ + L + C(jfc)[2 + |*] + C(^) 2 [2 + |] j. 


2 4 


Mh(k) = I + M 2 logMjy + ^C(m, 


(3.7) 


r- - 2 

M a (k) = ^ + i + M 2 logM|^ + y + C(i) 


a * HI 


For a = - y we have 


LJk) = k 2 l h (k), L a (k) = V~L a {k), M h (k) = ~rM h (k), M a (k) = —k ~M a (k). 


T 2- - 


- 2 , 


- 2 , 


Note: In the asymptotic expansions of the modified Bessel functions, we suppressed the dependence of y 
on the reduced frequency (it). The dependence is such that the higher the value of k, the smaller the 
value of M will need to be. However, for most aeroelastic calculations, we are only interested in the first 
few modes where \k\ is < 1 . In fact we have 


y = M 


Xb 

t/p 2 


l Xb 

P 2a ~ 






(3.8) 


4. COMPARISON WITH PREVIOUS RESULTS 


The definitions of the lift and moment coefficients in (3.6) are similar to those adopted by Possio, 
Dietze, and Timman. Timman’s lift and moment coefficients are defined for rotation about the midchord, 
a = 0, hence his L a . These moment coefficients are different from equation (3.6) and Dietze’s values. 
For M = 0.5, figures 3-6 compare the aerodynamic coefficients in equation (3.6) with their 
corresponding values found in Dietze and Timman. We also included the Prandtl-Glauert transformation 
of (3.6), i.e. all coefficients are divided by P . As k — > 0 , the Prandtl-Glauert version will reach the 
correct steady state value. Upon close inspection we found that Dietze s L ^ and L a were missing — k 
and -k 2 /2 respectively. When these errors are taken into account, our results show good agreement 
with Dietze’s values. The differences between Dietze’s values and ours are all within 10 percent. With 
the adjustment of the rotation axis, the differences between Timman’s values and ours are also within 
10 percent, except at high reduced frequencies. At these high reduced frequencies Fettis found that 
Timman’s values are erroneous. 16 As for the moment coefficients, again Dietze left out some terms there 
as well. When these errors are taken into account the difference between Dietze’s numbers and ours are 
again within 10 percent (not shown). 
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For the case where M = 0.7, figures 7-10 compare results from equation (3.7) with those that 
belong to Fettis 2 and Blanch 11 whenever they follow Smilg and Wassermann notations. Blanch 
employed Reissner’s method to calculate the aerodynamic coefficients, whereas Fettis approximated 
Possio’s equation in obtaining his results. It is uncanny that the values are within a few hundredths of 
each other. Our numbers compare reasonably well with theirs for Im[Ly,] and Re[L a ] shown in 
figures 8 and 9 respectively. The differences are close to 15 to 20 percent. However, for Re[L^] and 
Im[L a ] , there is nearly a 60-percent difference between our results and theirs. For this particular Mach 
number Fettis’ values agreed with Dietze’s when he took into account Dietze’s missing factors 
mentioned earlier. There were some difficulties in the convergence of Dietze’s iteration scheme for such 
a high Mach number. Fettis bypassed these difficulties by approximating Possio’s integral equation with 
a polynomial. He expanded the kernel function in a power series of k after he had removed the 
incompressible part. On the other hand, Blanch followed Reissner s method of expansions in terms of the 
Mathieu functions. Hence, their approximations are no longer valid for large Mach numbers and large 
reduced frequencies. Neither Fettis nor Blanch indicated how many terms should be taken in their 
expansions or what the limiting values are for the Mach number or the reduced frequency in their 
approximations. Our results are also approximations, but they are analytic rather than numerical 
approximations. If our approximations are good enough for Im[L^] and Re [ L a ] , they should be good 
enough for R e[L h ] and Im[L a ] as well. The differences are even greater in the moment coefficients. 
These differences are 60-percent and beyond (not shown). Just as with lift coefficients, we simply cannot 
conclude which approximation is better. From (3.7) it follows that 

M h (k) + i a (k ) = (l + i)l A (*), 

which is a formula derived by Fettis in reference 15. This relation is exact and was derived analytically 
based solely on the kinematics of the airfoil and the existence and uniqueness of the solution to the Possio 
equation. Overall, for M < 0.7 and k < 1 , our approximations are in good agreement with previous 
tabulated values. 
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Figure 7. Lift coefficient due to plunge motion as a function of reduced frequency ( k ) with M = 0.7 
(real part). 
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Figure 8. Lift coefficient due to plunge motion as a function of reduced frequency (k) with M = 0.7 
(imaginary part). 



Figure 9. Lift coefficient due to pitch motion as a function of reduced frequency ( k ) with M = 0.7 
(real part). 
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Figure 10. Lift coefficient due to pitch motion as a function of reduced frequency ( k ) with M = 0.7 
(imaginary part). 


5. APPLICATION TO FLUTTER ANALYSIS 


A main motivation for calculating the lift and moment in unsteady aerodynamics is for use in 
aeroelastic stability analysis — flutter speed calculations. In this section, we calculate the “aeroelastic 
modes” (see reference 13, p. 550) of a wing with a large aspect ratio and demonstrate the simplicity and 
the ease of utility of the explicit lift and moment formulas found in section 3. We will also see the need 
for the Laplace transforms of lift and moment as opposed to the “oscillatory motion” values when 
calculating these “aeroelastic modes.” 

We illustrate this with the continuum model of Goland 17 for bending-torsion flutter for a uniform 
cantilever wing. The flutter speed has been calculated in references 17 and 18 for the incompressible 
case. Here we shall use our formulas to calculate the aerodynamic loading in the subsonic compressible 
case. 


The dynamic equation for the wing is given as 

+ s + E/ fo_<y> = 

dr 2 ’ dt 2 dy* 

l d*g(v.O + 5 d 2 My,<) _ _ 

y . 2 y ^ 2 2 


ar 


ar dy' 

0<t, 0 <y<l. 


-L(y,t), 


(5.1) 
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with boundary conditions of. 


h(0,t) 


dhj 0,Q 

dy 


a(0 ,t) = 0, and 


_ a 3 ft(/,r) _ dajU) 

dy 2 dy 3 dy 


= 0 . 


Goland used the lift and moment of an infinitely spanned wing, but in actuality due to their dependence 
on h and a , the lift and moments in (5.1) do depend on y. For a wing with a sufficiently large aspect ratio, 
one can argue that the lift and moment can be written as the sum of the lift and moment of the infinitely 
spanned wing with a correction term that goes to zero as y — » 00 - 19 

Therefore, we proceed by taking the Laplace transform of equation (5.1) and denote the derivative 
with respect to y by superprime. We then have 


mX 2 h(y,X) + S y X 2 a(y,X) + EIh""(y,X) = -L(y,X), 
I^ 2 a(y,X) + S y X 1 h(y,'k)-GJa"(y,X) = M(y,X), 
0 < y < / 


with boundary conditions: 


h( 0,X) = h( 0,X) = a(0,X) = h"(l,X) = h”’(Uk) = a’(/,X) = 0. 


These equations yield the “aeroelastic” modes. The main point here is that we need L(y, X), M{ y, X) 
given by (3.6), (X is not purely imaginary). We may then proceed to calculate the real part of the 
aeroelastic mode corresponding to the first torsion mode, following reference 18, where the only 
difference is that the aerodynamics is restricted to the incompressible case, as in Goland. Figures 1 1 
and 12 plot the real part of X as a function of U, for the compressible as well as the incompressible case, 
at two different elevations. The flutter speed differs little from the incompressible case prior to flutter, 
even though the damping is higher initially. Tables 1 and 2 summarize the wing parameters and the 
flutter analysis. All of the reduced frequencies and Mach numbers at flutter condition are well within the 
range of validity of our approximation. 
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Figure 1 1 . Re [A] as a function of air speed: Goland case at sea level. 



Figure 12. Re [A.] as a function of air speed: Goland case at 20,000 ft. 





Table 1 . Parameters of the wing: Goland case. 


Parameters 

Values 

m 

0.746 slug/ft 

s y 

0.447 slug/ft/ft 


1.943 slug/ft 2 /ft 

El 

23.6 x 10 6 lb ft 

GJ 

2.39 x 10 6 lb ft 

a 

-1/3 

l 

20 ft 

b 

3 ft 


Table 2. Flutter analy sis of the wing: Goland case. 

Sea level 20 K ft above sea level 

Uf (ft/sec) k f M l f U f (ft/sec) kf My 

Incompressible 447 0.47 00 574 0.36 00 

Compressible 446 0.46 0.40 574 0.35 0.55 


6. CONCLUDING REMARKS 


For application to aeroelastic analysis in the subsonic compressible regime, the current reliance on 
tabulated calculations of aerodynamic influence coefficients is cumbersome. In this paper, closed-form 
formulas of lift and moment coefficients are obtained based on the recent solution of the Possio equation 
by Balakrishnan. 

For a Mach number less than 0.7 and reduced frequencies less than one, the closed-form formulas of 
lift and moment coefficients derived in this paper are within 10 percent of previous numerical results. 
As the Mach number gets closer to 0.7 and beyond, all methods of approximation fail. However, for such 
a high Mach number the aerodynamic flow is close to the transonic regime. 

Numerical calculations are shown to be consistent with extant work. The ease of use of these 
formulas is demonstrated in the flutter analysis illustrated. No longer is there a need for “table look up” of 
aerodynamic coefficients in an iterative process of root finding. 
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